REAL(8) FUNCTION func_gN2(x0)
! The values for gN and lump-sum transfers are preset as global variable, debt is optimally chosen

USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

REAL(8)::v2(na,nb2),c1(na,nb2),uc1(na,nb2)
REAL(8)::wmax,p0max,h0max,gn0max,ct0max,m0max,tau0max,ww0max,f0max,b_next_global
REAL(8)::value1,value2,value3
INTEGER  i_b, i_a, i_default_global,now_time
INTEGER::index_b,new_index_b,number_of_threads,i_b_base,try
INTEGER(4)::hours2,mins2,zero_h
INTEGER(4)::IRIGHT,ILEFT
REAL(8)::a_initial,b_initial,udef,u_value,indicator_external,gg
REAL(8)::stats(6),xlb(6),xub(6),penalty,time,ct_matrix_new(nb, na)
EXTERNAL udef
INTEGER(4), PARAMETER::nb3=501
REAL(8)::end_time_loop,start_time_loop      
REAL(8)::objective_u_fe,objective_u_unemp
EXTERNAL objective_u_fe,objective_u_unemp
REAL(8)::objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN_wbar,objective_u_fe_wbargN
EXTERNAL objective_u_fe_gN,objective_u_unemp_gN,objective_fun_gN_wbar,objective_u_fe_wbargN
REAL(8)::bgrid_base(nb3),h_baseX(na,nb3),gn_baseX(na,nb3),ct_baseX(na,nb3),bp_baseX(na,nb3),pn_baseX(na,nb3)    
REAL(8)::w_base(na,nb),pn_base(na,nb),ct_base(na,nb),aux0
REAL(8),DIMENSION(nb,na)::v_matrix_new,v0_matrix_new,vaut_matrix_new
REAL(8)::p0,term1,cnval,ctval 
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL,epsilongrid(ncdf), dev_q_nodef

index_bfeas2(:) = 0
iterfun = iterfun+1

call compute_bgrid

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)


!Initialize terminal node with default as risk-free debt environment
q_matrix_nodef = payoff/ (lambda + r)
v0_matrix   = 0d0
vaut_matrix = 0d0

do i_a = 1,na         
    FDATA(1:nb) = v0_matrix(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
	DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA(1:nb), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, 1:nb), coeff_matrix(i_a, :,1:nb))
ENDDO

do i_a = 1,na
    FDATA = q_matrix_nodef(:,i_a)
	ILEFT  = 0
	IRIGHT = 0
    DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
    break_matrix_q(i_a, :) = BREAK_GRID
    coeff_matrix_q(i_a, :,:) = CSCOEF
end do


open (214, FILE='output_flexwagesX\gh.txt')
open (215, FILE='output_flexwagesX\ggn.txt')
open (216, FILE='output_flexwagesX\gibp.txt')

open (314, FILE='output_flexwagesX\ghaut.txt')
open (315, FILE='output_flexwagesX\ggnaut.txt')
open (316, FILE='output_flexwagesX\bgrid2.txt')
open (317, FILE='output_flexwagesX\gp.txt')
open (318, FILE='output_flexwagesX\gct.txt')


DO i_b = 1,nb3
    READ(316,*) bgrid_base(i_b)    
ENDDO

CLOSE(316)

DO i_a = 1,na
    zero_h = 0
    DO i_b = 1,nb3
        READ(214,*) h_baseX(i_a,i_b)
        READ(215,*) gn_baseX(i_a,i_b)
        READ(317,*) pn_baseX(i_a,i_b)    
        READ(318,*) ct_baseX(i_a,i_b)    
 
        READ(216,'(I7)') i_b_base
        bp_baseX(i_a,i_b) = bgrid_base(i_b_base)    
        
        IF (h_baseX(i_a,i_b).EQ.0d0) THEN
            zero_h=i_b
        ENDIF
    ENDDO

    IF (zero_h.EQ.nb3) THEN
        PRINT*,'Problem: policy h not defined for any b given i_a=',i_a
        PAUSE
        STOP
    ENDIF

    zero_h = zero_h+1

    IF (zero_h.GT.1) THEN
        DO i_b = 1,zero_h-1
        gn_baseX(i_a,i_b) = gn_baseX(i_a,zero_h) !LINEAR_INT(bgrid_base(zero_h:101),gn_baseX(i_a,zero_h:101),bgrid_base(i_b))
        h_baseX(i_a,i_b)  = h_baseX(i_a,zero_h)  !LINEAR_INT(bgrid_base(zero_h:101),h_baseX(i_a,zero_h:101),bgrid_base(i_b))
            bp_baseX(i_a,i_b) = bp_baseX(i_a,zero_h)  !LINEAR_INT(bgrid_base(zero_h:101),h_baseX(i_a,zero_h:101),bgrid_base(i_b))
            pn_baseX(i_a,i_b) = pn_baseX(i_a,zero_h)  !LINEAR_INT(bgrid_base(zero_h:101),h_baseX(i_a,zero_h:101),bgrid_base(i_b))
            ct_baseX(i_a,i_b) = ct_baseX(i_a,zero_h)  !LINEAR_INT(bgrid_base(zero_h:101),h_baseX(i_a,zero_h:101),bgrid_base(i_b))
        ENDDO
    ENDIF
    READ(314,*) haut_base(i_a)
    READ(315,*) gnaut_base(i_a)
ENDDO

CLOSE(214)
CLOSE(215)
CLOSE(216)
CLOSE(314)
CLOSE(315)
CLOSE(317)
CLOSE(318)

     
DO i_a = 1,na
DO i_b = 1,nb
    gn_base(i_a,i_b) = LINEAR_INT(bgrid_base,gn_baseX(i_a,:),bgrid(i_b))
    h_base(i_a,i_b)  = LINEAR_INT(bgrid_base,h_baseX(i_a,:),bgrid(i_b))
    pn_base(i_a,i_b) = LINEAR_INT(bgrid_base,pn_baseX(i_a,:),bgrid(i_b))
    ct_base(i_a,i_b) = LINEAR_INT(bgrid_base,ct_baseX(i_a,:),bgrid(i_b))
    w_base(i_a,i_b) = alpha*pn_base(i_a,i_b)*h_base(i_a,i_b)**(alpha-1d0)
ENDDO
ENDDO

IF (psi1.LT.500d0) THEN
    open (230, FILE='output_flexwagesX\gqbase2.txt')
    open (231, FILE='output_flexwagesX\bgrid2base2.txt')
    DO i_a = 1,na
    DO i_b = 1,nb3
        READ(230,*) q1base(i_a,i_b)
    ENDDO
    ENDDO
    DO i_b = 1,nb3
        READ(231,*) bgridbase(i_b)
    ENDDO

    CLOSE(230)
    CLOSE(231)
ENDIF
    
indicator_external = 1 !FROM EXTERNAL FILE
if (indicator_external < 0.5) then
      !INITIAL BOUNDS ON GRID FOR b
    
        open (10, FILE='graphs\v.txt',STATUS='replace')
        open (11, FILE='graphs\default.txt',STATUS='replace')
        open (12, FILE='graphs\x.txt',STATUS='replace')
        open (13, FILE='graphs\gn.txt',STATUS='replace')
        open (14, FILE='graphs\b_next.txt',STATUS='replace')
        open (16, FILE='graphs\ctt.txt',STATUS='replace')

        do i_b = 1,nb 
            do i_a = 1,na 

                a_initial = agrid(i_a)
                b_initial = bgrid(i_b)
                try=1
                i_default_global = 2 !defaults  

                IF (choice_constant_gn.EQ.1)  constant_gn = gnaut_base(i_a)
                
                IF (psi1.LE.500d0) THEN
                    value1 = objective_u_fe(0d0,a_initial,b_initial,i_default_global)
                    value2 = objective_u_unemp(0d0,a_initial,b_initial,i_default_global)                             
                    
                    u_value = MAX(value1,value2)                 
                    vaut_matrix_new(i_b, i_a) = u_value-udef(a_initial)                    
                 ELSE
                    vaut_matrix_new(i_b, i_a) = -1d10                    
                 ENDIF
                  
                i_default_global = 1 !repays                 
                b_next_global=(1d0-lambda)*b_initial
                
                IF (choice_constant_gn.EQ.1)  constant_gn = gn_base(i_a,i_b)

                b_next_global=bmin
                
                value1 = objective_u_fe_gN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                value2 = objective_u_unemp_gN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                value3 = objective_u_fe_wbargN(b_next_global,constant_gn,a_initial,b_initial,i_default_global)
                
               IF (value1.GE.MAX(value2,value3)) THEN
                   x_matrix_new(i_b, i_a)=x_global1
                   gn_matrix_new(i_b, i_a)=gn_global1
                   ct_matrix_new(i_b, i_a)=ct_global1
               ELSEIF (value2.GE.MAX(value1,value3)) THEN
                   x_matrix_new(i_b, i_a)=x_global2
                   gn_matrix_new(i_b, i_a)=gn_global2
                   ct_matrix_new(i_b, i_a)=ct_global2
               ELSE
                   x_matrix_new(i_b, i_a) =x_global3
                   gn_matrix_new(i_b, i_a)=gn_global3
                   ct_matrix_new(i_b, i_a)=ct_global3
               ENDIF               
               
                u_value = MAX(value1,MAX(value2,value3))
                v0_matrix_new(i_b, i_a) = u_value  
                               
                v_matrix_new(i_b, i_a) = MAX(vaut_matrix_new(i_b, i_a), v0_matrix_new(i_b,i_a))

                 if (vaut_matrix_new(i_b, i_a) <= v0_matrix_new(i_b, i_a)) then
                   default_decision(i_b, i_a) = 1
                    b_next_matrix(i_b, i_a) = b_next_global
                 else
                   default_decision(i_b, i_a) = 2
                    b_next_matrix(i_b, i_a) = 0d0
                 end if              

          
                 WRITE(10, '(F18.10, X, F18.10, X, F15.10)') MAX(-100d0,v_matrix_new(i_b, i_a)),MAX(-100d0,v0_matrix_new(i_b, i_a)),MAX(-100d0,vaut_matrix_new(i_b,i_a))
                 WRITE(11, '(F18.10)') defaultgrid(default_decision(i_b, i_a))
                 WRITE(12, '(F18.10)') x_matrix_new(i_b, i_a)
                 WRITE(13, '(F18.10)') gn_matrix_new(i_b, i_a)
                 WRITE(14, '(F15.11)') b_next_matrix(i_b, i_a)
                 WRITE(16, '(F18.10)') ct_matrix_new(i_b, i_a)
                 
                 IF (psi1.LT.500d0) THEN        !default vs. risk-free debt economy         
                     q_matrix(i_b, i_a) = 0d0                 
                     q_matrix_nodef(i_b, i_a) = 0d0 
                 ELSE
                     q_matrix(i_b, i_a) = payoff/ (lambda + r)
                     q_matrix_nodef(i_b, i_a) = payoff/ (lambda + r)
                 ENDIF
                 
            end do
      end do
  
         CLOSE(10)
         CLOSE(11)
         CLOSE(12)
         CLOSE(13)
         CLOSE(14)
         CLOSE(16)

        v_matrix    = v_matrix_new
        v0_matrix   = v0_matrix_new
        vaut_matrix = vaut_matrix_new

                
        PRINT*, 'No Initial Guess Uploaded'
        PRINT*,' '
else !READ DATA FROM EXTERNAL FILES

open (10, FILE='graphs\v.txt')
open (16, FILE='graphs\q_paid.txt')
open (17, FILE='graphs\b_next.txt')
open (23, FILE='graphs\bounds.txt')

      READ(23, '(F15.11, X, F15.11)') bmin, bmax

      do i_b = 1,nb
         do i_a = 1,na

             READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrix(i_b, i_a), &
                     v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
             READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
             READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
             if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
                 default_decision(i_b, i_a) =2
             else
                 default_decision(i_b, i_a) =1
             end if             
         end do
      end do
CLOSE(10)
CLOSE(16)
CLOSE(17)
CLOSE(23)

PRINT*, 'Saved Initial Guess Uploaded'

end if


!call iterate_gN
!call pol_fcns_gN
!CALL INITIAL
!CALL fiscal_multi
!call MC_SIM(stats)

func_gN2 = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0


call system_clock(count=now_time)      !Option B
time  = (now_time-start_time)/rate_time/MAX(1,iterfun-1)
hours2 = INT(FLOOR(time/3600))
mins2  = FLOOR((time/3600-(hours2))*60d0)
PRINT*,'Iter = ',iterfun, 'Time per iteration',hours2,'hs',mins2,'mins'


PRINT*, 'MODEL SOLVED AND MC STATS COMPUTED'
RETURN

END FUNCTION func_gN2
!******************************************************************************************
REAL(8) FUNCTION objective_u_fe_gN2(x,gg,transf_star,a_initial2,b_initial2,i_default_global2)
!wage rigidity constraint is slack

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2,transf_star
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::bpval,wval,tauval,ct_global
REAL(8)::p0,m0,h0,q,b
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval
REAL(8)::aux_lhs,balt,aux0alt,minc,maxc,maxc2,maxc3
REAL(8)::term1,taxbar,bound1,bound2,new_value,old_value
REAL(8)::objective_u_feold,gn_global1old,x_global1old,ct_global1old 
REAL(8)::q2,aux02,p02,gnval2 ,temp1,temp2
INTEGER,PARAMETER::num=200
REAL(8)::ctemp(num),tolval,transf


INTEGER(4)::i,jj,im,iax,feas,error,niter2,ftype,try,index_opt,i_num
LOGICAL::flag00
EXTERNAL q_fun

objective_u_feold = -1d10
tolval = 1d-8

b0 = b_initial2
a0 = DEXP(a_initial2)

b = x

q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0           !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval= a0+aux0
ct_global = ctval   
IF (ctval.LE.0d0) THEN
    objective_u_fe_gN2= huge_n +1d3*(0d0-ctval)             
    if (choice_glob.EQ.1) THEN
        gn_global1 = 0d0
        x_global1  = 0d0
        ct_global1 = 0d0
    ENDIF
    RETURN
ENDIF

gnval = gg
cnval = 1d0-gnval
ctval = a0+aux0

p0   = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
wval = alpha*p0
tauval = p0*gnval-aux0

taxbar = taxrate*(a0+in*p0)
transf = taxbar-tauval

flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(gnval.GT.1d0).OR.(ctval.LE.0d0)) 

IF (.NOT. flag00) THEN
    pres  = p0
    hres  = 1d0
    ctres = ctval
    taures= tauval
    wwres = wval
    gnres = gnval
    feas  = 1d0
ELSE
    ctres   = 0d0    
    feas    = 0d0
    taures  = tauval
    hres    = 0d0
    pres    = p0
    gnres   = gnval
    wwres   = wval
ENDIF

!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(wwres.LT.wbar)*(wbar-wwres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty +1d6*((taures.GT.taxbar)*(taures-taxbar))
IF (taxrate.LE.1d0) penalty = penalty -1d9*((transf-transf_star)**2d0)
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres  = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres  = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

cres = cres - (cnres.LE.1d-6)*(-1d6)
cres = cres - (ctres.LE.1d-6)*(-1d6)

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cres,1d-8))
ELSE 
    utilc = MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma) 
ENDIF

IF (sigmag.EQ.1d0) THEN 
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n +penalty 

if (choice_glob.EQ.1) THEN
    gn_global1 = gnres
    x_global1  = 1d0
    ct_global1 = ctres
ENDIF

objective_u_fe_gN2 = wres

end FUNCTION objective_u_fe_gN2
!***************************************************************************
REAL(8) FUNCTION objective_u_unemp_gN2(x,gg,transf_star,a_initial2,b_initial2,i_default_global2)
!wage rigidity constraint binds

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2,transf_star
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=200
REAL(8)::htemp(num),aahat,aa1    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minh,maxh,temp1,taxbar,netbc,bound1,bound2,new_value,old_value,tolval
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt
REAL(8)::objective_u_unemp_gNold
INTEGER(4)::i_num_star,i_num_star2,index_optmax,root
real(8)::sign_new,sign_old,sign_two,ctemp2,difc,old_valuemax,transf
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10
tolval=1d-8
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x

q  = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))
aux_lhs   = aux0+taxrate*(a0+in*wbar/alpha)                          !always smaller than cT
!aux_lhs < 0 follows from  pn*gn = aux0+tax < 0 
!                          alpha*p*h**alpha=wbar*h < wbar

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF( (choice_fe.EQ.1).OR. (choice_nowbar.EQ.1)) THEN
    objective_u_unemp_gN2= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  (aux_lhs.LE.0d0) THEN
    balt = b+1d-6
    q    = q_fun(balt,a_initial2) 
    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_unemp_gN2= huge_n+1d3*(aux0alt-aux0)          
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

!case 1: h<1 and nonbinding tax cap
try = 1
gnval = gg

ftype = 199
minh  = (gnval+1d-4)**(1d0/alpha)
!Alternative: minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval           !b8
maxh  = 1d0

old_value = AFunction(ftype,a0,gnval,ct_global,minh)
old_valuemax = 10d+33
i_num_star=0

index_opt = 0
index_optmax = 0
difc = (maxh-minh)/(REAL(num)-1d0)
       
140 DO i_num=1+i_num_star,num
    ftype = 199
    htemp(i_num) = minh+(maxh-minh)*(REAL(i_num)-1d0)/(REAL(num)-1d0)     
    new_value = AFunction(ftype,a0,gnval,ct_global,htemp(i_num))
    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
        
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 282
    ENDIF
    old_value = new_value

ENDDO       

282 IF (index_opt.EQ.0) THEN
    bound1=htemp(MAX(index_optmax-1,1))
    bound2=htemp(MIN(index_optmax+1,num))
ELSE
    bound1=htemp(MAX(index_opt-1,1))
    bound2=htemp(MIN(index_opt,num))        
ENDIF

hval=BrentRoots(ftype,a0,gnval,ct_global,bound1,bound2,1d-10,1000,fval,niter2,error)  

IF ((error.EQ.1)) THEN 
    hval = 0d0
    IF (i_num_star.LT.num)  go to 140
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: More Iterations!'
    PAUSE
    STOP
ENDIF

cnval = hval**alpha-gnval
cnval = MAX(1d-6,cnval)
p0    = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)

1053 tauval = p0*gnval-aux0   
taxbar = taxrate*(a0+in*p0*hval**alpha)

transf=taxbar-tauval

flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0.0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty+1d6*((taures.GT.taxbar)*(taures-taxbar))
IF (taxrate.LE.1d0) penalty = penalty-1d6*((transf-transf_star)**2d0)
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.hres**alpha)*(gnres-hres**alpha)+(gnres.LT.0d0)*(0d0-gnres))
penalty = penalty +1d3*((hres.GT.1d0)*(hres-1d0)+(hres.LT.0d0)*(0d0-hres))
   
cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_unemp_gN2 = wres
if (choice_glob.EQ.1) THEN
    gn_global2 = gnres
    x_global2 = hres
    ct_global2 = ctres
ENDIF

IF (root.EQ.0) THEN
    x_global2old  = x_global2
    gn_global2old = gn_global2
    ct_global2old = ct_global2
    
    objective_u_unemp_gNold = objective_u_unemp_gN2
    root=1
ELSEIF (root.GT.0) THEN
    IF (objective_u_unemp_gN2.GT.objective_u_unemp_gNold) THEN
        objective_u_unemp_gNold= objective_u_unemp_gN2
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF
    root=root+1
ENDIF
    
IF (i_num_star.LT.num)     GO TO 140

IF (objective_u_unemp_gNold.GT.objective_u_unemp_gN2) THEN
    objective_u_unemp_gN2= objective_u_unemp_gNold
        
    if (choice_glob.EQ.1) THEN
        gn_global2 = gn_global2old    
        x_global2  = x_global2old 
        ct_global2 = ct_global2old
    ENDIF
ENDIF
    

RETURN

end function objective_u_unemp_gN2
!***************************************************************************
REAL(8) FUNCTION objective_u_fe_wbargN2(x,gg,transf_star,a_initial2,b_initial2,i_default_global2)
!wage rigidity constraint binds and h=1

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,gg,a_initial2,b_initial2,transf_star
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=200
REAL(8)::htemp(num),aahat,aa1,tolval    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minh,maxh,temp1,taxbar,netbc,bound1,bound2,new_value,old_value,transf
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10
tolval = 1d-8

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x

q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))
aux_lhs   = aux0+taxrate*(a0+in*wbar/alpha)                          !always smaller than cT
!aux_lhs < 0 follows from  pn*gn = aux0+tax < 0 
!                          alpha*p*h**alpha=wbar*h < wbar

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF (choice_nowbar.EQ.1) THEN
    objective_u_fe_wbargN2= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  (aux_lhs.LE.0d0) THEN
    balt = b+1d-6
    q = q_fun(balt,a_initial2)

    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_fe_wbargN2= huge_n+1d3*(aux0alt-aux0)            
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

gnval = gg
hval  = 1d0
p0    = wbar/alpha
cnval = 1d0-gnval
    
tauval = p0*gnval-aux0   
taxbar = taxrate*(a0+in*p0*hval**alpha)
transf=taxbar-tauval

IF (DABS(p0-(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)).GT.1d-6) p0=0d0

flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0.0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
IF (taxrate.LE.1d0) penalty = penalty+1d6*((taures.GT.taxbar+tolval)*(taures-taxbar))
IF (taxrate.LE.1d0) penalty = penalty-1d6*((transf-transf_star)**2d0)
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))
penalty=penalty+ 1d3*(wbar/alpha.GT.1d-6+(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))*(wbar/alpha - 1d-6-(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))
penalty=penalty+ 1d3*(wbar/alpha+1d-6.LT.(1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu))*((1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)-wbar/alpha-1d-6)

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))

IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_fe_wbargN2 = wres
if (choice_glob.EQ.1) THEN
    gn_global3 = gnres
    x_global3 = hres
    ct_global3 = ctres
ENDIF  

end function objective_u_fe_wbargN2
!******************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun_gN2(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV
REAL(8)::acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun
REAL(8)::gg,acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_unemp_gN2
REAL(8)::exp_v_next_aut,acum3,objective_u_fe_gN2
EXTERNAL DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_unemp_gN2,objective_u_fe_gN2

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa 
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

value1 = objective_u_fe_gN2(bb,constant_gn,transf_global,a_initial2,b_initial2,i_default_global2)
value2 = objective_u_unemp_gN2(bb,constant_gn,transf_global,a_initial2,b_initial2,i_default_global2)

u_value = MAX(value1,value2)
objective_fun_gN2 = -(u_value - ucost) - beta * exp_v_next

if (choice_glob.EQ.1) THEN
    IF (value1.GE.value2) THEN
        reg_global = 1
        x_global    = x_global1
        gn_global    = gn_global1
        ctrad_global = ct_global1  
    ELSE
        reg_global = 2
        x_global     = x_global2
        gn_global    = gn_global2
        ctrad_global = ct_global2  
    ENDIF    
ENDIF

end FUNCTION objective_fun_gN2
!******************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun_gN_wbar2(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV

REAL(8)::acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun
REAL(8)::gg,acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_fe_wbargN2
REAL(8)::exp_v_next_aut,acum3
EXTERNAL DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_fe_wbargN2

REAL(8)::value1,value2 ,u_value,ucost,udef

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa 
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)
u_value = objective_u_fe_wbargN2(bb,constant_gn,transf_global,a_initial2,b_initial2,i_default_global2)

objective_fun_gN_wbar2 = -(u_value - ucost) - beta * exp_v_next

end FUNCTION objective_fun_gN_wbar2
!********************************************************************************************
subroutine optimize_gN2(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt,index_opt2, index_vector(1), num,num2
parameter (num=10000) 
DOUBLE PRECISION :: b_next, v_value, objective_fun_gN2, new_value,  old_value,old_value2, b, vector(nb),&
                    b_next_grid(num),b_next_grid1(num), STEP, BOUND, XACC, b_next_guess,index_max
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun_gN2, DUVMIF
REAL(8)::b_next2,v_value2,bmax_aux

REAL(8)::objective_u_unemp_gN2,objective_u_fe_gN2,value3,value1,value2
EXTERNAL::objective_u_unemp_gN2,objective_u_fe_gN2

num2=num

b_next_inf = bmin
b_next_sup = bmax
index_max=0
 
old_value = 10d+33
do i=1,num2
   b_next_grid1(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun_gN2(b_next_grid1(i),a_initial2,b_initial2,i_default_global2)

   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid1(i)
      old_value = new_value
      index_max =1
  end if
   if ((DABS(old_value-new_value) <= 1d-6)) then
      index_opt2 = i
      old_value2 = new_value
   end if
end do
num2=num

b_next_inf = bp_global-0.025d0
b_next_sup = bp_global+0.025d0

303 do i=1,num2

   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun_gN2(b_next_grid(i),a_initial2,b_initial2,i_default_global2)

   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
      index_max=2
  end if
   if ((DABS(old_value-new_value) <= 1d-6)) then
      index_opt2 = i
      old_value2 = new_value
   end if
end do
num2=num

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSEIF (((DABS(old_value2 - old_value)) < 1d-6).AND.(index_opt2.NE.index_opt)) THEN   
    b_next  = b_next_grid(index_opt2)
    IF (index_max.EQ.1)     b_next  = b_next_grid1(index_opt2)
    v_value = -objective_fun_gN2(b_next,a_initial2,b_initial2,i_default_global2)
ELSE        
    IF (index_max.EQ.1)     STEP = (b_next_grid1(2) - b_next_grid1(1))*0.5
    call DUVMIF(objective_fun22, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun22(b_next)  
end if


CALL optimize_wbar2(b_next2,v_value2,a_initial2,b_initial2,i_default_global2)
value3=v_value2


IF (v_value2.GT.v_value) THEN
    v_value=v_value2
    b_next=b_next2
    
    reg_global = 2
    x_global     = x_global3
    gn_global    = gn_global3
    ctrad_global = ct_global3  
ENDIF

CONTAINS

REAL(8) FUNCTION objective_fun22(x)

REAL(8),INTENT(IN)::x

objective_fun22 = objective_fun_gN2(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
end subroutine
!**********************************************************************************
subroutine optimize_wbar2(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt,index_opt2, index_vector(1), num, num_tirar
parameter (num=300, num_tirar = 500)  !(num=300, num_tirar = 500)
DOUBLE PRECISION :: b_next, v_value, objective_fun_gN_wbar2, new_value, old_value,old_value2, b, vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun_gN_wbar2, DUVMIF

b_next_inf = bmin
b_next_sup = bmax

old_value = 10d+33
do i=1,num
   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun_gN_wbar2(b_next_grid(i),a_initial2,b_initial2,i_default_global2)
   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
   if (DABS(old_value-new_value) <= 1d-6) then
      index_opt2 = i
      old_value2 = new_value
  end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value
ELSEIF ((DABS(old_value2 - old_value)) < 1d-6) THEN       
    b_next  = b_next_grid(index_opt2)
    v_value = -objective_fun_gN_wbar2(b_next,a_initial2,b_initial2,i_default_global2)
ELSE        
    call DUVMIF(objective_fun22, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
    v_value = -objective_fun22(b_next)      
end if


CONTAINS

REAL(8) FUNCTION objective_fun22(x)

REAL(8),INTENT(IN)::x

objective_fun22 = objective_fun_gN_wbar2(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
end subroutine
!**********************************************************************************